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The change in contact angles due to the injection of low salinity water or any other wettability 
altering agent in an oil-rich porous medium is modeled by a network model of disordered pores 
transporting two immiscible fluids. We introduce a dynamic wettability altering mechanism, where 
the time dependent wetting property of each pore is determined by the cumulative flow of water 
through it. Simulations are performed to reach steady-state for different possible alterations in the 
wetting angle (6). We find that deviation from oil-wet conditions re-mobilizes the stuck clusters and 
increases the oil fractional flow. However, the rate of increase in the fractional flow depends strongly 
on 9 and as 9 — > 90°, a critical angle, the system shows critical slowing down which is characterized 
by two dynamic critical exponents. 

PACS numbers: 47.56.+r,47.61.Jd 


I. INTRODUCTION 

The world’s primary energy demand is predicted to in¬ 
crease by one-third between 2011 and 2035, where 82% 
of it comes from fossil fuels [1], In this scenario, the fact 
that some 20 to 60 percent of the oil remains unrecovered 
in a reservoir after the production is declared unprof¬ 
itable, is a challenge of increasing importance [2]. The 
main reason for this loss is the formation of oil clusters 
trapped in water and held in place by capillary forces, 
which in turn are controlled by the wetting properties 
of the reservoir fluids with respect to the matrix rock. 
The production from complex oil reserves that today are 
considered immobile or too slow compared to the cost is 
therefore an important area of research. In this context, 
the role of formation wettability is a focus area within 
the field of Enhanced Oil Recovery (EOR) [3] . 

Different reservoir rocks have widely different wetting 
characteristics [4] . Wettability may vary at the pore level 
from strongly oil wet through intermediate wetting to 
strongly water wet. Carbonate reservoirs contain more 
than half of the world’s conventional oil reserves, but 
the oil recovery factor is very low compared to sandstone 
reservoirs [5]. This is due to the complex structure, for¬ 
mation heterogeneity and more chemically active wetta¬ 
bility characteristics of the carbonate reservoirs, which 
leads to uncertainty in the fluid flow and oil recovery [6] . 
Sandstone is strongly water wet before oil migrates from 
a source rock into the reservoir. When oil enters a pore, 
it displaces the water which leaves behind a water film 
sandwiched between the oil and rock surface. This hap¬ 
pens as a result of balancing van der Waals and electric 
double layer forces, capillary pressure and grain curva¬ 
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ture [7]. A permanent wettability alteration is then be¬ 
lieved to take place by adsorption of asphaltenes from the 
crude oil to the rock, and leads to high but slow recovery 
through continuous oil films [8, 9]. As the oil saturation 
drops, these films can become discontinuous, leaving im¬ 
mobile oil clusters held in place by capillary forces. 

Wettability is therefore an important petrophysical 
property which plays a key role in the fluid transport 
properties of both conventional and unconventional reser¬ 
voirs [10], and there is great potential to improve the 
oil recovery efficiency by altering the wetting properties. 
Main factors which can alter the pore wettability are; 
lowering the salinity [11], adding water-soluble surfac¬ 
tants [12, 13] or adding oil-soluble organic acids or bases 
[14]. Increasing the reservoir temperature also increases 
water-wetness [4, 15]. There are some correlations with 
the wetting behavior to the electrostatic forces between 
the rock and oil surfaces [16], but there is no consensus 
on the dominating microscopic mechanism behind the 
wettability alteration. It is known from laboratory ex¬ 
periments and field tests that a drift from strongly oil- 
wet to water-wet or intermediate-wet conditions can sig¬ 
nificantly improve the oil recovery efficiency [14]. The 
amount of change in the wetting angle is a key factor here 
[17, 18] which not only decides the increase in oil flow but 
also the speed of the process. An improper change in the 
wetting angle can also make the recovery very slow and 
not profitable. 

Given there is a certain change in the wetting angle due 
to a brine, the next important factor is the flow pathways 
in the matrix rock which transports the oil and brine. 
One cannot expect any change in the wetting angle of 
a pore if there is no flow of the brine through it. The 
flow pathways depend on several different factors: the 
porous network itself, oil saturation, capillary number 
and also on the present wettability conditions. A change 
in the wettability will cause a perturbation in the flow 
distribution of the system. This will in turn again affect 
the wettability change through the altered flow pathways, 
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causing further changes in the flow distribution. The 
dynamics of wettability alterations is therefore controlled 
by a strongly correlated process. 

There are some studies of wettability alterations in 
two-phase flow by equilibrium-based network models [19] 
for capillary dominated regimes where viscous forces are 
negligible. Moreover, investigation of the time-scale of 
dynamics lacks attention in such models. In this arti¬ 
cle, we present a detailed study of wettability alterations 
in two-phase flow considering a network model of dis¬ 
ordered pores transporting two immiscible fluids where 
a dynamic wettability alteration mechanism, correlated 
with the flow-pathways, is implemented. We will focus 
on the transport properties due to the change in the wet¬ 
tability as well as on the time scale of the dynamics. 

We study in the following the effect of wettability 
changes on immiscible two-phase flow based on a net¬ 
work model [21-23]. In section II, we present the model 
and how we adapt it to incorporate the dynamic wet¬ 
tability changes. In Section III, we present our results. 
Initially, we let the two phases settle into a steady state 
where the averages of the macroscopic flow parameters 
no longer evolve. At some point, we then introduce the 
wettability altering agent, so that it starts changing the 
wetting angle. The wetting angle alteration depends on 
the cumulative volume of the wettability altering fluid 
that has flowed past a given pore. This induces transient 
behavior in the macroscopic flow properties and we mea¬ 
sure the time it takes to settle back into a new steady 
state. We find that there is a critical point at a wetting 
angle of 90 degrees and we measure its dynamical critical 
exponents; the exponents are different whether one ap¬ 
proaches the critical point from smaller or larger angles. 
In Section IV we summarize and conclude. 


II. MODEL 

We model the porous medium by a network of tubes 
(or links) oriented at 45 degrees relative to the overall 
flow direction. The links contain volumes contributed 
from both the pore and the throat, which then intersect 
at volume-less vertices (or nodes). Any disorder can be 
introduced in the model by a proper random number 
distribution for the radius r of each link, and we choose a 
uniform distribution in the range [0.1Z, 0.4Z] here, where 
l is the length of each tube. The network transports two 
immiscible fluids (we name them as oil and water), one of 
which is more wetting than the other with respect to the 
pore surface. The pores are assumed to be in between 
particles, and the pore shape is thus approximated to be 
hour-glass shaped, which introduces capillary effects in 
the system. The model is illustrated in figure 1. 

Due to the hour glass shape of the pore, the capillary 
pressure at a menisci separating the two fluids is not 
constant, and depends on the position x of the menisci 
inside the pore. The capillary pressure p c (x ) at position 
x inside the ith pore is then calculated from a modified 



FIG. 1. (a) Illustration of the pore network model, con¬ 

structed by links (tubes) with random radii connected to each 
other via nodes at the intersection of the dashed lines. One 
single link in the network is highlighted by gray which is again 
shown in (b) filled with two fluids. The wetting and non¬ 
wetting fluids are colored by white and gray respectively, pi 
and p 2 are the pressures at the two ends of the link and q is the 
local flow-rate. There are three menisci inside and capillary 
pressure difference at a menisci is indicated as p c (x). 


form of the Young Laplace equation [20, 21], 
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where 7 is the interfacial tension between the fluids and 
is the wetting angle for that pore. As an interface 
moves in time, p c {x) changes. The capillary pressure is 
zero at the two ends (x = 0 and Z) and it is maximum 
at the narrowest part of the pore. It makes the model 
closer to the dynamics of drainage dominated flow, where 
the him flow can be neglected. When there are multi¬ 
ple menisci in a pore, the total capillary pressure inside 
the ith pore is obtained from the vector sum {Y^iPc(x)) 
over all the menisci in that pore. The flow is driven by 
maintaining a constant total flow rate Q throughout the 
network, which introduces a global pressure drop. The 
instantaneous local flow rate qi inside the ith link be¬ 
tween two nodes with pressures p\ and P 2 follows the 
Washburn equation of capillary flow [24], 
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where A pi = P2 — pi- ki = rf /8 is the permeabil¬ 
ity of cylindrical tubes. Any other cross-sectional shape 
will only lead to an additional overall geometrical factor. 

— PoSi + Pvj( 1 — Si), is the volume averaged viscosity 
of the two phases inside the link, which is a function of 
the oil saturation Si in that link. Here p, Q and /z w are the 
viscosities of oil and water respectively. 

The flow equations for the tube network are solved 
using a conjugate gradient method [25]. These are the 
Kirchhoff equations balancing the flow, where the net 
fluid flux through every node should be zero at each time 
step, combined with the constitutive equation relating 
flux and pressure drop across each tube. The system of 
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equations is then integrated in time using an explicit Eu¬ 
ler scheme with a discrete time step and all the menisci 
positions are changed accordingly. Inside the zth tube, 
all menisci move with a speed determined by qt- Physical 
processes like bubble snap-off and coalescence are intro¬ 
duced in the model. Due to these, bubbles can be formed 
or merged, which changes the number of menisci inside 
a link with time. When a meniscus reaches the end of a 
tube, new menisci are formed in the neighboring tubes 
depending upon the flow rates. In each link, a maximum 
number of menisci is allowed to form. When this number 
is exceeded, the two nearest menisci are merged, keeping 
the volume of each fluid conserved. In our simulations, 
we considered a maximum of four menisci inside one pore 
which can be tuned depending on the experimental obser¬ 
vations. Further details regarding the menisci movement 
can be found in Reference [22]. 

The simulations are started with an initial random dis¬ 
tribution of two fluids in a pure oil-wet network. Bi- 
periodic boundary conditions are implemented in the sys¬ 
tem, which effectively makes flow on a torus surface. The 
flow can therefore go on for infinite time, keeping the sat¬ 
uration constant and the system eventually reaches to a 
steady state. In the steady state, both drainage and im¬ 
bibition take place simultaneously and fluid clusters are 
created, merged and broken into small clusters. One can 
consider this as the secondary recovery stage. Once the 
system reaches the steady-state in a oil-wet network, the 
dynamic wettability alterations are implemented, which 
may be considered as the tertiary recovery stage or EOR. 
In the following we discuss this in detail. 


Dynamic wettability alteration 


We now introduce a dynamic wettability alteration 
mechanism to simulate any wetting angle change, decided 
by the oil-brine-rock combination and the distribution of 
the flow channels in the system. In a previous study 
[23] , a simplified static wettability alteration mechanism 
was studied, where the alteration probability was con¬ 
sidered equal for all pores without any correlation with 
the flow of brine inside a pore. However, for wettabil¬ 
ity alterations to occur, the wettability altering agent 
(e.g. low-salinity water or surfactant) needs to be in con¬ 
tact with the pore walls. Thus, the wettability alteration 
should follow the fluid flow pathways and any change in 
the wetting angle inside a pore should depend on the 
cumulative volumetric flow of brine through that pore. 
This claim is rather trivial, as one can not expect any 
wettability change in a pore if the altering agent is not 
present. This means that if a certain pore is flooded by 
large amounts of brine, the wetting angle should change 
more in that pore than the one which had very little wa¬ 
ter flooded through. This is implemented in the model 
by measuring the cumulative volumetric flux V) (t) in each 



FIG. 2. Change in the wetting angle for a link given by Gi(t) 
as a function of the cumulative volume of water Vi (t) through 
that link. 


individual pore with time t, 

t 

Vi(t) = ^2qi(t){l-s(t))At, (3) 

t=to 


where to is the time when the injection of low salinity wa¬ 
ter is initiated, At is the time interval between two sim¬ 
ulation steps and (1 — s(t)) is the water saturation. V)(t) 
is then used to change the wetting angle for each tube 
continuously, updated at every time step after t = to- 
The wetting angle of the zth pore can change contin¬ 
uously from 180° to 0° as V)(f) changes from 0 to oo. 
Correspondingly, the cosfii term in Eq. (1) will change 
from —1 to 1 continuously. This continuous change of 
the wetting angle with the variation of Vi ( t ) is modelled 
by a function Gi ( t ) given by, 


GAt) = — tan 

7T 


< 7(^-1 

V th 


( 4 ) 


which replaces the cos term in Eq. (1). The pre-factor 
2 /7T is a normalization constant to set the range of the 
function, and the parameter C adjusts the slope during 
the transition between oil wet and water wet. This leads 
to the change in the wetting angle as a function of V % ( t ) 
as shown in Fig. 2. 

As our model does not include film flow, the wetting 
angle should not reach either 0 or 180 degrees. Rather, a 
starting point of « 165° was used by choosing C = 20 
in our simulations as shown in Fig. 2. Next, as a larger 
pore will need more brine to be flooded in order to have a 
similar change in the wetting angle than a smaller pore, a 
threshold value V t th is introduced, which is proportional 
to the volume of that pore, 

Vt = rpnjl. (5) 

At Vi (t) = R/ h , the wetting angle reaches to 90° in that 
pore and p c (x) essentially becomes zero. Here p is a pro¬ 
portionality constant which decides how many pore vol¬ 
umes of water is needed to reach I/[(f) = V) th for the zth 
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FIG. 3. (Color online) Change in the oil-fractional-flow ( F) 
during the simulation as a function of the fluid volume passed 
through any cross-section of the network expressed in terms 
of the number of pore-volumes (IV). The network is of size 
40 x 40 links with oil saturation S = 0.3 and Ca = 10 -1 here. 
The initial part of the plot (N < 200) shows the change in 
fractional flow in an oil-wet system where it approaches to a 
steady state with F « 0.235. The wettability alteration starts 
at IV = 200, where results for simulations with different values 
of 77 are plotted in different colors. The system evolves to a 
new steady-state with a higher average fractional flow, which 
is same for any value of 77 . 77 only has an effect on the rate of 
change in F, which is shown in the inset. A higher value of rj 
results in a longer time to reach the new steady-state. 


pore. This parameter can possibly be adjusted against 
future experimental results, but is considered as a tuning 
parameter in this study. The expression for the capillary 
pressure at a menisci from Eq. 1 then takes the form, 


Pc(t ) 


2 ^Gj(t) 
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1 — cos 



( 6 ) 


The maximum amount of wetting angle that can be 
changed depends on the combination of brine, oil and 
rock properties [17, 18]. We therefore set a cut-off 9 in 
the wetting angle change, such that any pore that has 
reached to a wetting angle = 9 , can not be changed 
further. The model thus includes all the essential ingre¬ 
dients of wettability alteration study - it is a time depen¬ 
dent model where the wettability alteration is correlated 
with the flow pathways of the brine, and can be used to 
study any oil-brine-rock combination decided by 9. 


III. RESULTS 

Simulations are started with a random initial distri¬ 
bution of oil and water in an oil-wet network, where 
9 ss 165° for all links. First, the oil-wet system is 
evolved to a steady state before any wettability alteration 
is started. This will allow us to compare the change in 
the steady-state fractional flow of oil ( F ) with a change 
in the wetting angle. The oil fractional flow (F) is de¬ 
fined as the ratio of the oil flow-rate (Q 0 n) to the total 



FIG. 4. Distribution of fluid bubbles (top row) and local flow 
rates (bottom row) in steady state in a network of 64 x 64 links 
with oil saturation S = 0.3 and Ca = 10~ 2 . The left column is 
for the steady state in a oil-wet network and the right column 
is the new steady state after wettability alteration takes place. 
In (a) and (b), the oil bubbles are drawn in black. In (c) and 
(d), the normalized local flow-rates Qi/^max are drawn in gray 
scale. 


flow-rate (Q) given by, F = Q 0 n/Q- The flow rate (Q) is 
kept constant throughout the simulation, which sets the 
capillary number Ca = /r e f jQ/ (jA), where A is the cross- 
sectional area of the network. A network of 40 x 40 links 
are considered, which is sufficient to be in the asymp¬ 
totic limit for the range of parameters used [22], An av¬ 
erage over 5 different realizations of the network has been 
taken for each simulation. As the simulation continues, 
both drainage and imbibition take place simultaneously 
due to bi-periodic boundary conditions and the system 
eventually evolves to a steady state, with a distribution 
of water and oil clusters in the system. In Fig. 3, F is 
plotted against the number of pore volumes passed (IV) 
through the network. As we run the system with con¬ 
stant flow-rate, N is directly proportional to the time t, 
N = tQ/v where v is the total volume of the network. 
The initial 200 pore volumes are for an oil-wet network, 
where it reaches to a steady state with F ~ 0.235. We 
then initiate the dynamic wettability alteration which re¬ 
sembles the flow of a wettability altering brine and F 
starts to drift. Here we run simulations for different val¬ 
ues of 77 , defined in Eq. 5, and the results are plotted 
in different colors. 9 = 0° in these simulations, which 
means any pore can change to pure water-wet depending 
upon the flow of brine through it. One can see that F 
approaches to a new steady-state with F ~ 0.308 due to 
the wettability alteration. Here different values of 77 make 
the simulation faster or slower to reach the same steady 
state, where a higher value of 77 makes the system slower 
to reach the steady state. The initialization of steady 
state measured in terms of the pore-volumes (N ss ) is de¬ 
fined as the instant when the average fractional flow stops 
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FIG. 5. Oil fractional-flow ( F) in the steady state at different 
oil saturation (S) in the initial oil-wet system (O) and after 
the wettability alteration with 9 C = 0° (□). Individual simu¬ 
lations are performed for all S values at two different capillary 
numbers (a) Ca = 10 -1 and (b) Ca = 10 -2 . F is higher after 
wettability alteration for the whole range of S. 


changing with time anymore and essentially stays within 
its fluctuation. In the inset of Fig. 3, N* s is plotted with 
r] and a simple linear dependency is observed. In order 
to save computational time, we therefore use 77 = 10 in 
all our following simulations. 

How the two fluids and the local flow-rates are dis¬ 
tributed in the network in the two steady-states before 
and after the wettability alteration are shown in Fig. 
4. The network size is 64 x 64 links here with an oil- 
saturation S = 0.3 and the capillary number Ca = 10 -2 . 
All the links are hour-glass shaped in the actual sim¬ 
ulation with disorder in radii, but shown as a regular 
network for simplicity in drawing. The upper row shows 
the distribution of oil bubbles drawn in black. The left 
column (a) shows the steady state in a oil-wet network 
and the right column (b) shows the steady-state after the 
wettability alteration is initiated with maximum possible 
wetting angle change 9 = 0°. A closer look in these bub¬ 
ble distributions shows more clustered oil bubbles in (a) 
than in (b) where they are more fragmented. A more in¬ 
teresting picture can be seen in the local flow-rate distri¬ 
bution in the bottom row, which shows a more dynamic 
scenario. The left (c) and right (d) figures are for the 
same time-steps before and after wettability alteration 
as in (a) and (b) respectively. Here the local flow-rates 
in each pore, normalized in between 0 to 1 , are shown in 
gray scale. Interestingly, in the oil-wet system (c), the 
flow is dominated in a few channels (black lines) where 
the flow-rates are orders of magnitude higher than the 
rest of the system. Other than those channels, the sys¬ 
tem has negligible flow, indicated by white patches which 
means the fluids are effectively stuck in all those areas. 
This situation happens when the difference in the satura¬ 
tion of the two fluids is large, where the phase with higher 
saturation (water here) tries to percolate in paths dom¬ 
inated by a single phase with less number of interfaces. 
This is not favorable in oil-recovery, as it leaves immobile 
fluid in the reservoir. In the flow distribution after the 
wettability alteration (d), the flow is more homogeneous 
and distributed over the whole system, indicating higher 
mobility of the fluids. However, one should remember 
that when the wettability alteration is started in a sys- 



FIG. 6. Proportionate change in the steady-state oil 
fractional-flow (A F/F) due to wettability alteration as a func¬ 
tion of maximum wetting angle 9 for (a) Ca = ICG 1 with 
S = 0.3 and (b) Ca = 10 -2 with S = 0.4. 


tern shown in (c), the wettability alteration starts tak¬ 
ing place only in those pores with active flow. But then 
it perturbs the flow-field and starts new flow paths and 
eventually the system drifts towards a more homogeneous 
flow with time, as shown in (d). 

We now present the results when the wetting angle of 
any pore can change all the way down to zero degree 
(9 = 0°). In Fig. 5 the steady-state oil fractional-flow 
in an initial oil-wet system ( F ) is compared with that in 
the steady-state after wettability alteration (F'). Results 
are plotted as a function of S for two different capillary 
numbers, (a) Ca = 10 _1 and (b) Ca = 10 -2 . The diago¬ 
nal line corresponds to F = S. A miscible fluid mixture 
would follow this line, and it is interesting to use this 
as a reference to compare how the data points lie above 
or below this line. If both the fluids flow equally in the 
system then F will be equal to S. But in Fig. 5, F stays 
below the F = S line for low oil saturation values in the 
oil-wet system, which means that the flow of oil is less 
than that of water and there are stuck region of oil. A 
significant increase in F can be observed here due to the 
wettability alteration for the full range of oil-saturation. 
Moreover, increase in F is higher for the lower capillary 
number, indicating that wettability alteration is very sig¬ 
nificant in the case of oil recovery, as Ca can go as low 
as 10 ~ 6 in the reservoir pores. Fractional flow also obeys 
the symmetry relation F'(S) = 1 — F( 1 — S) [23] which 
implies that, if the wetting angle of any pore is allowed to 
change all way down to zero degree (0 = 0 °), the system 
will eventually become pure water-wet with time. 

As noted earlier, the maximum change in wetting angle 
for a system, depends on the properties of the reservoir 
rock, crude oil and brine, and also on the temperature. 
Existing wettability alteration procedures generally turns 
the oil-wet system into intermediate wet, rather than to 
pure water-wet. Some examples of the change in the wet¬ 
ting angle for different rock materials and brine can be 
found in [17, 18]. In our simulation this is taken care of 
by the parameter 9 , which decides the maximum change 
in the wetting angle for any pore. One should re¬ 
member that, we are not forcibly changing the wetting 
angles 1 9, to 9 , rather the change in is decided inde¬ 
pendently for individual pores by the amount of brine 
passed through it (by Eqs. 3 to 6 ), and there is a max- 
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FIG. 7. Change in F during the simulation for different 9 
values. The wettability alteration started at N = 400 pore- 
volumes. The initialization of the steady-states for different 
9 values are marked by crosses on the respective plots. 



FIG. 8. Variation of steady-state initialization time (r) as 
a function of 9 for (a) Ca = 10 -1 with S = 0.3 and (b) 
Ca = 10 -2 with S = 0.4. r diverges rapidly as 9 —> 9 C where 
9 C = 90°. 

imum allowed change in any $i. As before, simulations 
are started with a pure oil-wet system to reach a steady 
state and then wettability alteration is started and simu¬ 
lation continues until the system reaches to a steady state 
again. Independent simulations have been performed for 
different values of 9. The proportionate change in the 
oil fractional-flow due this wettability change from the 
oil-wet system, A F/F = (F' — F)/F is measured for 
different simulations and plotted in Fig. 6 against 9. 
There are a few things to notice. First, as one can im¬ 
mediately see, fractional flow increases with the decrease 
of oil-wetness, 9 —> 0°. The maximum increase in F is 
higher for lower Ca, about 86% for Ca = 10~ 2 and about 
32% for Ca = 10 -1 . This is because the change in wet¬ 
ting angle affects the capillary pressures at the menisci, 
so the change in F is larger when the capillary forces are 
higher. Secondly, the major change in F happens in the 
intermediate wetting regime, upto 9 ss 60°, and then it 
becomes almost flat afterwards. Moreover, this increase 
in F is more rapid for lower value of Ca. All these facts 
points towards an optimal range of wetting angle change 
to increase the oil flow. This is an important observation 
for practical reasons, as it is not necessary to change the 
wetting angle further. Thirdly, there is a discontinuity in 
the curve exactly at 9 = 90°, as we will discuss later. 

Increase in the oil fractional-flow with the increase in 
the water-wetness may seem to be obvious and reported 


by many experiments and held tests. But, the most im¬ 
portant concern for the oil industry is the rate of increase, 
or the time required to achieve a significant increase in 
the oil production. If the increment in oil flow is very 
slow compared to the cost of the process, then the oil 
recovery is declared as not profitable and the reservoir 
may be considered to be abandoned. As per our knowl¬ 
edge, there are very few systematic studies reported in 
the literature predicting the time scale to change the oil 
how due to the wettability change by two-phase how of 
brine and oil in a porous media. We observe that, due 
to the correlations between the how paths and the wet¬ 
ting angle change, the time scale of the process varies 
dramatically with 9. This is illustrated in Fig. 7, where 
F is plotted as a function of pore volumes ( N ) of huid 
passed through the system. The initial 400 pore volumes 
are for an oil-wet system and then results of few differ¬ 
ent simulations with 9 = 85, 88, 89, 90, 91, 92 and 94 
degrees are plotted. Interestingly, the rate at which the 
system reaches a new steady-state, varies signihcantly 
depending on the value of 9. For example, after the wet¬ 
tability alteration is started, it needs to how less than 100 
pore volumes to reach the new steady state for 9 = 94° 
whereas more than 300 pore volumes are needed to reach 
a steady state for 9 = 91°. Therefore, even if the hnal 
steady-state fractional how is higher for 9 = 91° than for 
9 = 94°, it might not be profitable to alter the wetting 
angles to 91° because of the slow increase in F. In gen¬ 
eral, the process becomes slower and slower as 9 —> 90° 
from both sides. Such kind of slow increase in oil recov¬ 
ery as 9 90° is also observed in experiments [17, 18]. 

This slowing down of the process is an combined effect 
of two factors. First, the fact that wettability only can 
change in the pores where there is how of brine and the 
second is the value of 9. All the pores were initially oil- 
wet (i?j « 165°) and when it reaches the steady state, 
the how hnds the high mobility pathways depending on 
the mobility factor of the pores and the capillary pres¬ 
sures at the menisci. When the wettability alteration 
is started, the wetting angles of the existing how path¬ 
ways start decreasing. As a result, capillary pressures at 
menisci in those channels first decreases as 90° and 
then it increase afterwards as —>• 0°. This creates a 

perturbation in the global pressure held and correspond¬ 
ingly viscous pressures start changing with time which 
changes the how held. However, capillary pressures at 
the zero-how regimes are now higher than the high-how 
regimes which makes it difficult to invade the zero-how 
regimes causing a slower change in the how held as 
approaches 90°. An interesting feature is observed ex¬ 
actly at 9 = 90°, where the average fractional how does 
not change at all after the wettability change. At exactly 
9 = 90°, capillary pressures in all the pores in the exist¬ 
ing how pathways essentially become zero, making them 
the lowest resistive channels with zero capillary barri¬ 
ers. As a result, the huids keep howing in the existing 
channels forever and the system stays in the same steady- 
state. The time taken to reach another new steady-state 
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FIG. 9. Plot of logr vs log|0 — 0 C \ for (a) Ca = 10 _1 with 
S = 0.3 and (b) Ca = 10~ 2 with S = 0.4. From the slopes, 
the value of the dynamic critical exponent a is obtained as 
a = 1.0±0.1 for 9 < 9 C (Q) and a = 1.2±0.1 for 9 > 6 C (□). 


is therefore infinite at 0 = 90° and it therefore is a critical 
point for the system. 

We now measure the steady-state initialization time t, 
defined as the moment when the average of the fractional 
flow stops changing with time and becomes horizontal 
with the x axis. This is shown In Fig. 7, where the 
initialization of steady states is marked by crosses on 
the respective plots. As the simulations are performed 
with constant Q, r is proportional to the fluid volume 
passed through the system and therefore we measure r in 
the units of N. r for different simulations with different 
values of the maximum wetting angle (0) is plotted in Fig. 

8 for (a) Ca = 10 -1 with S = 0.3 and (b) Ca = 10~ 2 
with S = 0.4. One can see that r diverges rapidly as 

9 approaches 9 C = 90° from both sides, 9 > 90° and 

< 90°. This divergence of the steady-state time r as 6 —> 
9 C indicates the critical slowing down of the dynamics, 
which is a characteristics of critical phenomena. The 
critical slowing down is the outcome of the divergence of 
correlations at the critical point and can be characterized 
by a dynamic critical exponent z defined as r ~ , where 

£ is the correlation length [26]. As 9 —> 9 Cl the correlation 
length £ diverges as \9—9 c \~ 1 ' and therefore r ~ \9—9 c \~ a , 
where a = zv. In Fig. 9, r is plotted versus \9—9 c \ in log- 
log scale which gives two different slopes for 9 > 9 C and 
9 < 9 C . We then find the value of the dynamic exponents 
a as a = 1.2 ± 0.1 for 9 > 9 C and a = 1.0 ± 0.1 for 
9 < 9 C . However, they are the same within error bar for 
different capillary numbers and saturations (Fig. 9). The 
value of the dynamic critical exponents depend on the 
underlying dynamics and on the model [27] . In this case, 


wettability alteration was started from an oil-wet system 
with di = 165° for all the pores. So for the simulations 
with 9 < 90°, the wetting angles cross the critical point 
(90°) when the capillary forces change directions. This 
might cause the system to mobilize the clusters somewhat 
faster than for 9 > 90° when the capillary forces does not 
change any direction. As a result, a becomes smaller for 
9 < 90° than for 9 > 90°. 

IV. CONCLUSIONS 

In this article we have presented a detailed computa¬ 
tional study of wettability alterations in two-phase flow 
in porous media, where the change in the wetting angle in 
a pore is controlled by the volumetric flow of the altering 
agent through it. When the wetting angles are allowed to 
alter towards water-wetness, the stuck oil clusters start to 
mobilize and oil-fractional flow increases. However, due 
to the correlations in the wetting angle change with the 
flow pathways, the time-scale of the dynamics strongly 
depends on the maximum allowed change in the wetting 
angle. We find that, as the final wetting angle is chosen 
closer to 90°, the system shows a critical slowing down 
in the dynamics. This critical slowing down is charac¬ 
terized by two dynamic critical exponents. The critical 
point we are dealing with is an equilibrium critical point 
as the system is in steady state. The dynamical criti¬ 
cal exponents measure how long it takes to go from one 
steady state to a new one. To our knowledge, this is the 
first example of there being different values for the expo¬ 
nents on either side of the critical point. Our findings are 
in agreement with experimental observations reported in 
literature, and are extremely important for application 
purposes like oil recovery, where the time scale of the 
process is a key issue. 
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